#reading in fits data file
from astropy.coordinates import SkyCoord
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.table import Table
import numpy as np
I’ve put the dataset on Canvas, but it’s also good to see what the MAST is and what it can do for you (hint, it’s amazing!).
Here are the directions for downloading the data from MAST (please read even if you get data off of Canvas):
Go to the Space Telescope archive (http://archive.stsci.edu/mast.html) and download GHRS ECH-B spectra of the star HZ 43. What kind of a star is HZ 43? (Hint, consult SIMBAD). The GHRS dataset contains a series of .fits files. c0f.fits contains the wavelength solutions, c1f.fits contains the calibrated fluxes, *c2f.fits contains the uncertainties. Consult http://www.stsci.edu/documents/dhb/pdf/GHRS.pdf for details. There are two HST GHRS ECH-B datasets for HZ 43. Choose the one that contains the Mg II absorption lines near 2796 and 2803 ̊A. This dataset appears to have 4 identical copies available on the archive - two from HST online and two from DADS.
#reading in fits files
#wavelength solutions
wavedata= fits.getdata('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c0f.fits')
print(wavedata)
headerwave= fits.getheader('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c0f.fits')
print(wavedata.shape)
#calibrated fluxes
calfluxdata= fits.getdata('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c1f.fits')
print(calfluxdata)
#uncertainities
uncertdata= fits.getdata('/Users/sandratredinnick/Desktop/ASTR5550/z2r50205t.c2f.fits')
print(uncertdata)
[[2792.38625404 2792.39387204 2792.40148979 ... 2807.30055332 2807.3078941 2807.31523656] [2792.38625404 2792.39387204 2792.40148979 ... 2807.30055332 2807.3078941 2807.31523656] [2792.38625404 2792.39387204 2792.40148979 ... 2807.30055332 2807.3078941 2807.31523656] ... [2791.96152474 2791.96914634 2791.97676771 ... 2806.88311864 2806.89046314 2806.89780933] [2791.96152474 2791.96914634 2791.97676771 ... 2806.88311864 2806.89046314 2806.89780933] [2791.96152474 2791.96914634 2791.97676771 ... 2806.88311864 2806.89046314 2806.89780933]] (16, 2000) [[2.2250016e-13 3.0200262e-13 6.3603928e-14 ... 1.7262567e-13 5.5018222e-14 2.1968450e-13] [1.4303592e-13 1.4305775e-13 1.4307957e-13 ... 3.3727434e-13 2.6671747e-13 2.4320300e-13] [3.7150437e-13 2.9208875e-13 5.3690247e-14 ... 1.8643812e-13 1.3939940e-13 2.8054525e-13] ... [1.3972841e-13 2.1924210e-13 1.3977200e-13 ... 2.6314199e-13 9.8672908e-14 1.9266676e-13] [1.3913449e-13 2.1864817e-13 2.1868071e-13 ... 1.9207281e-13 1.6858249e-13 1.9208588e-13] [3.7853298e-13 2.9909487e-13 1.4013291e-13 ... 9.9022244e-14 1.6951636e-13 1.9301975e-13]] [[1.37627336e-13 1.58939851e-13 7.94806658e-14 ... 6.65300876e-14 4.07424779e-14 7.43876261e-14] [1.12372263e-13 1.12387462e-13 1.12402648e-13 ... 9.11000740e-14 8.14849558e-14 7.80183955e-14] [1.77676152e-13 1.58939851e-13 7.94806658e-14 ... 7.05658134e-14 6.22351698e-14 8.48149404e-14] ... [1.12403434e-13 1.37683986e-13 1.12433649e-13 ... 8.13928935e-14 5.25405186e-14 7.04926704e-14] [1.12403434e-13 1.37683986e-13 1.37702539e-13 ... 7.04883200e-14 6.64590724e-14 7.04926704e-14] [1.77725429e-13 1.58983774e-13 1.12433649e-13 ... 5.25388957e-14 6.64590724e-14 7.04926704e-14]]
According to SIMBAD, star HX 43 is a White Drawf!
Look through the GHRS.pdf document to understand why the 16 subexposures in the dataset have four different sets of wavelength solutions (look in Section 35.6). Follow the recommended procedure and align the subexposures by cross-correlating the fluxes (here, don’t just use numpy’s or scipy’s correlate functions!) against the first exposure obtained, then shifting by the number of wavelength bins required to achieve maximum cross-correlation. You may want to co-add each set of 4 before cross-correlating to improve s/n. You should get more than 10 pixel shifts, and you have to think about what the correlation does, to do it properly. Eq 3.102 of Ivezic should help connect the integral definition in class to the discrete function.
#plot data first to see what im working with
#plot first row of wavelengths vs flux
##16 rows so it seems like those are my 16 sub exporesure
plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[0],calfluxdata[0])
plt.xlabel("Wavelength")
plt.ylabel("Flux")
Text(0, 0.5, 'Flux')
plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[1],calfluxdata[1])
plt.xlabel("Wavelength")
plt.ylabel("Flux")
Text(0, 0.5, 'Flux')
plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[2],calfluxdata[2])
plt.xlabel("Wavelength")
plt.ylabel("Flux")
Text(0, 0.5, 'Flux')
plt.figure(figsize=(20,10), dpi=400)
plt.plot(wavedata[3],calfluxdata[3])
plt.xlabel("Wavelength")
plt.ylabel("Flux")
Text(0, 0.5, 'Flux')
#finding the 4 groups of the same wavelenth (so i can co-add spectra)
## co-add: mean of the flux arrays 0-3,4-7,8-11,12-15
print(wavedata[0]==wavedata[3])
print(wavedata[4]==wavedata[7])
print(wavedata[8]==wavedata[11])
print(wavedata[12]==wavedata[15])
[ True True True ... True True True] [ True True True ... True True True] [ True True True ... True True True] [ True True True ... True True True]
#wavedata[0]==wavedata[2]
#make array of 4 sets to co add spectra
array1= np.array(calfluxdata[0:4])
array2= np.array(calfluxdata[4:8])
array3=np.array(calfluxdata[8:12])
array4= np.array(calfluxdata[12:16])
#co adding spectra
setone= np.mean(array1,axis=0)
settwo= np.mean(array2,axis=0)
setthree= np.mean(array3,axis=0)
setfour= np.mean(array4,axis=0)
To create the cross validation function in python, I will use the following equation which was pulled from wikipedia:
The idea of this function is to sweep the second flux array over the first flux array from right to left. The overlapping area between the two flux arrays is the cross correlation. Np.roll() inverts arrays. First I need to invert the first flux array then multiply each element of the first flux array by each element of the second flux array (sweeping across). Then need to sum over that mulitplication to get cross correlation.
#n= index
#f = flux array one
#g = flux array two
#sum over all values in m
#shift by m
####GHRS.pdf document said:The best way to align the subexposures
#is by cross-correlating them
#against the first exposure obtained
#create cross correlation function because cant use np or scipy correlate funct.
N=100
def crosscorr(spectra_one,spectra_two):
crosscorr=np.zeros(N)
for i in range (N):
x= np.roll(spectra_two,-i) #np.roll inverts the second array flux
y=x*spectra_one #multiply each element above by element of other flux array
crosscorr[i]=sum(y) #sum over the multiplication
return crosscorr
To create the cross correlation function I worked with Anna.
#cross correlate the fluxes against the first exposure obtained
ccfone= crosscorr(setone,settwo)
ccftwo=crosscorr(setone,setthree)
ccfthree= crosscorr(setone,setfour)
#shift by the number of wavelength bins
#required to achieve maximum cross-correlation
#find the shift, number of wavelength binds required to achieve max correlation
shiftone = np.argmax(ccfone)
shifttwo = np.argmax(ccftwo)
shiftthree = np.argmax(ccfthree)
print(shiftone)
print(shifttwo)
print(shiftthree)
17 33 52
# allign/shifting subexposures
allignone= np.roll(settwo,-shiftone)
alligntwo= np.roll(setthree,-shifttwo)
allignthree= np.roll(setfour,-shiftthree)
The pixel shifts or number of wavelength bins required to ahcieve maximum cross-correlation are (17, 33, and 52).
Co-add (i.e., find the mean of) all aligned subexposures, using uniform weighting since presumably each of the sub-exposures has the same integration time. Calculate the uncertainty in your co-added spectrum using the given uncertainties for the individual spectra. Turn in a plot of the co-added spectrum, with the uncertainty values overplotted in some fashion (either along the bottom of the plot or plt.errorbar). Make a plot of the s/n as a function of wavelength. Make sure to label axes with the correct units.
#co adding all alligned subexposures
arraycoadd=np.array([setone,allignone,alligntwo,allignthree])
coadd2= np.mean(arraycoadd,axis=0) #this is my master flux array
#calculating uncertainities in co adding
#need to use propogation of errors --> boy am I am getting good at prop of errors!
unccoaddone = ((uncertdata[0])**2+(uncertdata[1])**2+(uncertdata[2])**2+(uncertdata[3])**2)/16
#print(unccoaddone)
unccoaddtwo = (uncertdata[4]**2+uncertdata[5]**2+uncertdata[6]**2+uncertdata[7]**2)/16
#print(unccoaddtwo)
unccoaddthree = (uncertdata[8]**2+uncertdata[9]**2+uncertdata[10]**2+uncertdata[11]**2)/16
#print(unccoaddthree)
unccoaddfour = (uncertdata[12]**2+uncertdata[13]**2+uncertdata[14]**2+uncertdata[15]**2)/16
#print(unccoaddfour)
Because I didn't want to latex my hand written work, see canvas dropbox submission for written work for this question.
# shifting
unccoaddshiftone = np.roll(unccoaddtwo, -shiftone)
unccoaddshifttwo = np.roll(unccoaddthree, -shifttwo)
unccoaddshiftthree = np.roll(unccoaddfour, -shiftthree)
#finalized calculated uncertainity values
sigma2_unccoaddshift= (unccoaddone+ unccoaddshiftone+ unccoaddshifttwo+unccoaddshiftthree)/16
sigma_unccoaddshift = np.sqrt(sigma2_unccoaddshift)
#plot co added spectrum with uncertainities
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavedata[0],coadd2,label="Spectrum")
plt.errorbar(wavedata[0],sigma_unccoaddshift,label="Uncertainity")
plt.fill_between(wavedata[0],coadd2-sigma_unccoaddshift,coadd2+sigma_unccoaddshift,color='black',label="Uncertainity")
plt.title("Co-added Spectrum With Associated Uncertainity Values")
plt.xlabel("Wavelength ($\AA$)")
plt.ylabel("Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)")
plt.legend()
<matplotlib.legend.Legend at 0x7fe278a01310>
#find s/n
#s = master flux array
#n= errors
#recall: coadd2 is the alligned co added fluxes
master_flux = coadd2
signoise= master_flux/sigma_unccoaddshift
print(signoise)
print(np.mean(signoise))
[ 9.59806 9.402988 9.870239 ... 12.560999 12.236406 13.652048] 12.339193
#plot s/n
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavedata[0],signoise)
plt.xlabel("Wavelength ($\AA$)")
plt.ylabel("Signal to Noise Ratio (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)")
plt.title("Spectrum of Signal to Noise Ratio as a Function of Wavelength")
#plt.axvline(2796.3, color = 'black', label='Mg II absorption line')
#plt.axvline(2803.48, color = 'black', label='Mg II absorption line')
Text(0.5, 1.0, 'Spectrum of Signal to Noise Ratio as a Function of Wavelength')
#need to find indicies of absoption lines!
#note the absoption lines are a range, not just
#a single horizontal line
np.argmin(signoise),1000+ np.argmin(signoise[1000:])
(515, 1479)
Is the s/n between the two absorption lines consistent with the standard deviation of the flux values? What is the minimum percent change relative to the continuum between the lines at which you could confidently detect an emission or absorption line?
# find the average continuum flux value between
#the two Mg II lines, and divided by the standard deviation in that region
#find the signal to nosie array between the abs lines
#so did 515+5, 1479-5 to get inbetween array
signoisebetween = signoise[520:1474]
#mean contiuum value beetween the 2 abs lines
meansignoisebetween = np.mean(signoisebetween)
print(meansignoisebetween)
#std vaue of between the 2 abs line contiuum
stdsignoisebetween= np.std(signoisebetween)
print(stdsignoisebetween)
12.324845 0.61785144
#find s/n between the two absorption lines for flux
#not using 515:1480 bc the abs line is a range of values not
#just a single point
fluxbetween= master_flux[520:1474]
#mean of the the between contiuum
meanfluxbetween= np.mean(fluxbetween)
print(meanfluxbetween)
# find the standard devition of the between continuum
stdfluxbetween= np.std(master_flux[520:1474])
print(stdfluxbetween)
#find the s/n for the contiuum between the 2 abs lines
print(meanfluxbetween/stdfluxbetween)
2.326488e-13 1.9952174e-14 11.660323
I am assuming that you want us to see if the s/n between the absorption lines is consistent with the mean +/- standard deviation of the flux values.
I found that the s/n between the absorption lines is consistent with the mean +/- the standard deviation of the flux values. I got the mean +/- the standard deviation of flux values to be 12.32 +/- 0.62. I got the s/n between the two absorption lines to be 11.66. Thus, yes the s/n between the two absorption lines is consistent with the mean +/- the standard deviation of the flux values.
percent = (abs((meanfluxbetween-3*stdfluxbetween-meanfluxbetween)/meanfluxbetween))*100
print(percent)
25.728274398227423
The minimum percent change relative to the continuum between the lines at which I could confidently detect an emisson or absorption line is approximately 25.73%.
Because I didn't want to latex my hand written work, see canvas dropbox submission for written work for this question.
Search for the interstellar absorption lines of Mg II (near 2796 and 2803 ̊A). What are the laboratory wavelengths of these lines? (See Morton, D. 1991, ApJS, 77, 119 for this and oscillator strengths (f values); also can use NIST database: https://physics.nist.gov/PhysRefData/ASD/lines_form.html and make sure to select the appropriate wavelength scale (air or vacuum).)
According to NIST database, the laboratory wavelengths in a vacuum scale of these absorption lines are 2796.352 ̊A and 2803.530 ̊A.
What is the s/n in the continuum on either side of these lines? What is the s/n in the absorption line? If the noise is dominated by Poisson noise from photon counts, what would you expect the relationship between fractional uncertainties in the continuum and the line to be? Check to see if your s/n ratio is consistent with a Poisson distribution.
#find s/n contintuum (array) on either side of the absoption lines
#goal: find s/n on left and right sides of absoption lines
#so need mean and standard deviation of left and right sides! (s/n)
#recall the absoption line is a range of values not a single point
#first--> left side
#find array that is the left side of asborption line
spectrumleft= coadd2[:510]
#finding mean of spectra to the left side of absoption line
spectrumleftmean = np.mean(spectrumleft)
#findind std of spectra to the left side of absoption line
spectrumleftstd = np.std(spectrumleft)
#print(spectrumleftmean)
#print(spectrumleftstd)
#s/n for left side of asbotion lines (0 to 514)
spectrumleftSN = spectrumleftmean/spectrumleftstd
print(spectrumleftSN)
#second --> right side
#find array that is the right side of asborption line
spectrumright= coadd2[1480+5:]
#finding mean of spectra to the right side of absoption line
spectrumrightmean = np.mean(spectrumright)
#findind std of spectra to the right side of absoption line
spectrumrightstd = np.std(spectrumright)
#print(spectrumrightmean)
#print(spectrumrightstd)
#s/n for right side of asbotion lines (1480 to end)
spectrumrightSN = spectrumrightmean/spectrumrightstd
print(spectrumrightSN)
12.267691 12.0922575
The s/n for the continuum on the left side of the left MgII absorption line (line near 2796 ̊A) is 12.23. The s/n for the continuum on the right side of the right absorption line (line near 2803 ̊A) is 10.84. And, from part 1.4 we saw that the s/n for the continum for the right side of the left MgII absorption line (line near 2796 ̊A) and the left side of the ight absorption line (line near 2803 ̊A) is 12.09.
#find the S/N in two absorption line
#first find the S/N in the left absoption line (line near 2796 ̊A)
#signoise = signal to noise of whole continum so just
#get s/n of the absoption line
leftabslineSN = signoise[515]
print(leftabslineSN)
#second find the s/n for the right asboption line (line near 2803 ̊A)
rightabslineSN = signoise[1479]
print(rightabslineSN)
5.1924796 7.073334
The s/n for the left asboption (near 2796 ̊A) is 5.19. The s/n for the right absoption line (near 2803 ̊A) is 7.07.
#fractional uncertainity stuff
#goal: assuming the noise is dominated by Poisson noise from photon counts
##need to compute fractional uncertainities in the continuum and the line
#use def of fractional uncertanity, and to look at the relationship between them
#take the ratio --> fractional uncertainities in continuum/fractional unct in line
#note that uncertainities for poisson stats are defined by sqrt(mu)!
#deciding to the continuum inbetween the absorption lines
#calculating the ratio for the line, chooshing the left absorption line
ratiocontline = np.sqrt(coadd2[515])/np.sqrt(np.mean(coadd2[520:1469]))
print(ratiocontline)
0.47227946
#find s/n ratio for the conintuum inbetween the 2 abs lines
SNratio = (np.std(coadd2[520:1469])*coadd2[515])/(np.mean(coadd2[520:1469])*sigma_unccoaddshift[515])
print(SNratio)
#now can compare!
0.44459233
I expect the relationship/ratio between the fractional uncertainties in the continuum between the absorption lines and the left absoprtion line (near 2796 ̊A) to be 0.472. The s/n ratio I got for the continuum inbetween the two absorption lines and using the left absorption line (near 2796 ̊A) was 0.445. The values are pretty close, so yes, my s/n ratio is consistent with a Poisson distribtuion.
Because I didn't want to latex the formulas needed for this question (assuming noise is dominated from a poisson distribtuion, consider relationship between fractional uncertainties in the continuum and the line), see canvas dropbox submission for written work for this question.
Fit a low order polynomial (even just a straight line with a slope is fine) to the continuum, masking out the absorption lines. What are the uncertainties in your fit parameters? Turn in a plot of the spectrum with your fit to the background overplotted. (This could be a good chance to write up that design matrix algorithm I’ve been talking about, but a canned routine is fine too.).
Coni helped me figure out how to mask out the absoption lines (i.e. delete the points in the array that were at the absoption lines). I tried using np.remove(), as well as add arrays to left, right, and inbetween the absoption lines (but arrays werent of the same shapes), and more, but getting lots of errors that I couldn't debug quickly enough. She sent me her code and I used it as you said I could over email. I understand completely how the code works:
#first step: mask out the absorption lines
spectrumleftnew= np.where(wavedata[0]<wavedata[0][500])[0]
spectrumbetweennew= np.where((wavedata[0]>wavedata[0][530]) & (wavedata[0]<wavedata[0][1454]))[0]
spectrumrightnew = np.where(wavedata[0] > wavedata[0][1494])[0]
indicies = [spectrumleftnew,spectrumbetweennew,spectrumrightnew]
indicies2 = np.concatenate(indicies)
meanspectrummask= coadd2[indicies2]
wavesmask= wavedata[0][(indicies2)]
#plot the spectrum or continuum of the two absoprtion
#lines masked out
#note that the line inbetween is Python artifically connecting
#the data points but the line doesnt have meaning
plt.plot(wavesmask, meanspectrummask, color='black')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.title('Spectrum With Asborption Lines Masked Out')
plt.show()
#doing a polynomial fit
#doing a degree one polynomial fit
#array for params and fit
parameters,fit = np.polyfit(wavesmask,meanspectrummask,deg=1,cov=True)
#creating an array for the fit to make a smooth curve
x = np.linspace(wavesmask[0],wavesmask[-1],1000)
#creating the line fit
##parameters[0]=slope
y = parameters[0]*x + parameters[1]
#plot masked spectrum with fit overplotted
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavesmask, meanspectrummask, color='black',label='Specturm')
plt.plot(x,y,label='Degree One Polynomial Fit')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.legend()
#plt.axhline(np.mean(meanspectrummask), c='green',label='Mean of spectrum')
plt.title('Spectrum With Asborption Lines Masked Out & Polynomial Fit')
plt.show()
#can fix code to be different, take sqrt inside
#uncertainities in fit paramters
sig2slopeone = fit[0][0]
sig2inttwo = fit[1][1]
#uncertainity in slope paramater
sigslopeone = np.sqrt(sig2slopeone)
#uncertainity in intercept paramater
siginttwo = np.sqrt(sig2inttwo)
print(sigslopeone)
print(siginttwo)
#print paramter estimates
#slope
print(parameters[0])
#intercept
print(parameters[1])
1.0115119626338277e-16 2.832114711126018e-13 2.8235668811163627e-16 -5.583116058824341e-13
I fit a degree one polynomial to the continuum which masked out the two absorption lines. Note that a degree one polynomial fit makes the model: y= slope*x+intercept. The parameter estimates with their corresponding uncertainities for this degree one polynomial fit are:
1. Intercept: -5.583116058824341e-13 +/- 2.832114711126018e-13
2. Slope: 2.8235668811163627e-16 +/- 1.0115119626338277e-16
Normalize the spectrum by dividing by your best fit continuum polynomial or line. The continuum should now be relatively flat, with a mean value of unity in the continuum. Turn in a plot of the normalized spectrum.
#creating best fit line for normalized specturm, abs lines
#masked out
bestfit = parameters[0]*wavesmask+parameters[1]
#normalizing the spectrum with abs lines masked out
normspec= meanspectrummask/bestfit
normspec
#sanity checking to make sure it is normalized.
#if it was normalized would get mean =1
meannorm = np.mean(normspec)
print(meannorm)
0.999999889775847
#plotting the normalized spectrum
plt.figure(dpi = 200, figsize = (8, 6))
plt.plot(wavesmask, normspec, color='black')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.title('Normalized Spectrum With Asborption Lines Masked Out')
#plt.axhline(meannorm, c='green',label='Mean of spectrum')
plt.show()
I confirm that my the continuum is relatively flat, with a mean value of unity in the continumum.
Fit each of the two absorption lines with a Gaussian. Report your best-fit values and uncertainties for the fit parameters, including the central wavelength, amplitude, and width. Are the single Gaussian fits acceptable? Turn in a plot of the normalized spectrum with the Gaussian fits overplotted.
#goal: do gaussian fit to both absorption lines
#step one: first need normalized spectrum with
#the asbroption lines in
###so I must renormalize with a best fit line that considers
#the whole spectrum instead of the specturm with the absroption
#lines removed
#same procedure as 2.4
#getting best fit line, using whole spectrum with
#absorption lines in
bestfit2 = parameters[0]*wavedata[0]+parameters[1]
#normalized spectrum array with abs lines in
normspec2= coadd2/bestfit2
#sanity check--> make sure its actually normalized by checking mean
#mean should be close to 1!
meannorm2 = np.mean(normspec2)
#plot Normalized Spectrum/continuum with Abs Lines in spectrum
plt.plot(wavedata[0], normspec2, color='black')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.title('Normalized Spectrum With Absoption Lines')
#plt.axhline(meannorm, c='green',label='Mean of spectrum')
plt.show()
#now I need to find two different arrays where the two different arrays
#are the first absoption line and second abs
#####similar process to 2.3, but flipped, want only the absoption lines
#note I am using the same code as 2.3 to get these two arrays (the code i used
#from 2.3 was Conis and I am applying it again here, see 2.3 for my explanation
#of exactly what the code is doing)
#first absorption line array
lineone= np.where((wavedata[0]>wavedata[0][505]) & (wavedata[0]<wavedata[0][525]))[0]
print(lineone)
#first absoption line flux array and wavelength array
##note that I multiplied the flux array by -1 because I need to invert the flux array
#to properly fit a gaussian
#fluxarray
fluxlineone= -normspec2[lineone]+1 #shift up by one to get on right scale
#wavelength array
waveslineone= wavedata[0][lineone]
[506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524]
#plot the absoption line continuum for the first line (line to the left of the spectrum)
#plt.scatter(waveslineone, fluxlineone, color='black')
#plt.xlabel('Wavelength')
#plt.ylabel('Flux')
#plt.title('Spectrum of Absorption Line One (Near 2796 $\AA$)')
#plt.show()
To fit the gaussian I copied code from:
https://stackoverflow.com/questions/59047395/fitting-gaussian-to-absorbtion-line-in-python
from scipy.optimize import curve_fit
#gaussian fit for the absoption lines
## I found this code from stack overflow! (see link directly above)
#parameters of fit are: central wavelength (mean),amplitude, and width (standard deviation)
#first define the gaussian function
def Gauss(x, amp, mean, sigma):
return amp * np.exp(-(x - mean)**2 / (2 * sigma**2))
#x= wavelength array of the first abs line
#amp= amplitude
#need starting values for the gaussian fit of abs line one!
#fit by eye
#worked together with coni to set the starting values
am= np.max(fluxlineone)-np.min(fluxlineone)
mean= waveslineone[np.argmin(fluxlineone)]
std=.09
#doing the fitting for absorption line 1
param,unc = curve_fit(Gauss, waveslineone, fluxlineone, p0=[np.max(fluxlineone)-np.min(fluxlineone),waveslineone[np.argmin(fluxlineone)],.09])
#want to create an array of wavelengths to make the gaussian fit smooth for abs line 1
x1= np.linspace(min(waveslineone),max(waveslineone),1000)
#plotting the asboption line one continuum/spectrum with fit
plt.plot(x1,-Gauss(x1,param[0],param[1],param[2])+1,c='red') #invert the gaussian fit
#function wont fit a inverse gaussian! sad! but working around it.
plt.scatter(waveslineone, -fluxlineone+1, color='black') #invert the fluxes
plt.title('Normalized Spectrum of Absorption Line One (Near 2796 $\AA$) With Gaussian Fit')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.show()
#best fit values and uncertanities for fit parameters for abs line one
#best fit values for amplitude, mean, standard deviation respectively
print(param)
#uncertainties associated with amplitude, mean, standard deviation respectively
#the diagnoials are the uncertanties associated with each paramater
print(unc)
[7.83477723e-01 2.79629393e+03 2.80868715e-02] [[ 6.62321913e-04 -1.59543516e-07 -1.62450183e-05] [-1.59543516e-07 1.13232367e-06 1.49759747e-08] [-1.62450183e-05 1.49759747e-08 1.17286212e-06]]
The diagoanals are the uncertainties associated with each parameter, respectively.
#second absorption line array
linetwo= np.where((wavedata[0]>wavedata[0][1469]) & (wavedata[0]<wavedata[0][1489]))[0]
print(linetwo)
#first absoption line flux array and wavelength array
##note that I multiplied the flux array by -1 because I need to invert the flux array
#to properly fit a gaussian
fluxlinetwo= -normspec2[linetwo]+1
waveslinetwo= wavedata[0][linetwo]
[1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488]
#plot the absoption line continuum for the second line (line to the left of the spectrum)
#plt.scatter(waveslinetwo, fluxlinetwo, color='black')
#plt.xlabel('Wavelength')
#plt.ylabel('Flux')
#plt.title('Spectrum of Absorption Line Two (Near 2803 $\AA$)')
#plt.show()
#need starting values for the gaussian fit for abs line 2 !
#fit by eye
#worked together with coni to set the starting values
am3= np.max(fluxlinetwo)-np.min(fluxlinetwo)
mean3= waveslinetwo[np.argmin(fluxlinetwo)]
std3=.02
#doing the fitting for absorption line 2
param2,unc2 = curve_fit(Gauss, waveslinetwo, fluxlinetwo, p0=[np.max(fluxlinetwo)-np.min(fluxlinetwo),mean3,std3])
#want to create an array of wavelengths to make the gaussian fit smooth for abs line 1
x2= np.linspace(min(waveslinetwo),max(waveslinetwo),1000)
#plotting the asboption line two continuum/spectrum with fit
plt.plot(x2,-Gauss(x2,param2[0],param2[1],param2[2])+1,c='red') #invert the gaussian fit
#function wont fit a inverse gaussian! sad! but working around it.
plt.scatter(waveslinetwo, -fluxlinetwo+1, color='black') #invert the fluxes
plt.title('Normalized Spectrum of Absorption Line Two (Near 2803 $\AA$) With Gaussian Fit')
plt.xlabel('Wavelength ($\AA$)')
plt.ylabel('Flux (erg cm$^{-2}$ s$^{-1}$ $\AA^{-1}$)')
plt.show()
I fit the two absorption lines with a Gaussian. The parameters used to fit the model were central wavelength (mean), amplitude, and width (standard deviation).
The parameter estimates for the second fit of the first absorption line (near 2796 ̊A) with their corresponding uncertainities are as follows:
1. Mean: 2.79629393e+03 +/- 0.0010641069824035552
2. Ampltiude: 7.83477723e-01 +/- 0.02573561565224349
3. Standard Deviation: 2.80868715e-02 +/- 0.0010829875899565978
The parameter estimates for the Gaussian fit of the second absorption line (near 2803 ̊A) with their corresponding uncertainities are as follows:
1. Mean: 2.80347993e+03 +/- 0.0018593463851579672
2. Amplitude: 6.19312544e-01 +/- 0.0400935179299597
3. Standard Deviation: 2.49001433e-02 +/- 0.0018673348494579111
Yes, the Gaussian fits of the absorption lines are acceptable. The plot of both of the two fits looks pretty good and the uncertainities are super small! I'm shocked at how well they fit the data.
Calculate the central velocity for each line, and the uncertainties in the central velocities. Use the convention that a positive velocity is away from the observer (careful!). Are the central velocities of these lines consistent with each other? Test the consistency by using a differenced χ2 consistency test, taking into account the uncertainties in the laboratory wavelengths of these lines.
#goal: calculate the central velocity for each line
#will need to use propogation of errors
#best fit line one
###amp,mean,std=[7.83477723e-01 2.79629393e+03 2.80868715e-02]
##uncetainiies: 6.62321913e-04, 1.13232367e-06, 1.17286212e-06]
#best fit line two
##amp,mean,std= [6.19312544e-01, 2.80347993e+03, 2.49001433e-02]
##uncertainities: 1.60749018e-03, 3.45716898e-06, 3.48693944e-06
#observed = parameter values found in 3.1
#central velcoties for abs line one
lambdarestone= 2796.352 #found in 2.1
light= 2.99e+8 #speed of light
centvelone= light*(param[1]-lambdarestone)/lambdarestone # in m/s
#errors for velecoties for abs lines one
##use prop of errors
uncvelone = (light/lambdarestone)*np.sqrt(unc[1][1])
#observed = parameter values found in 3.1
#central velcoties for abs line two
lambdaresttwo= 2803.530 #found in 2.1
centveltwo= light*(param2[1]-lambdaresttwo)/lambdaresttwo
#errors for velecoties for abs lines 2
##use prop of errors
uncveltwo = (light/lambdaresttwo)*np.sqrt(unc2[1][1])
print(centvelone)
print(uncvelone)
print(centveltwo)
print(uncveltwo)
-6209.12674621733 113.77966292471372 -5340.074961801732 198.30163008107883
The central velocity with its associated uncertainities for the first absorption line (near 2796 ̊A) is:
The central velocities for the two absorption lines are relatively differently, thus the central velocities for the two absorption lines are inconsistent with each other.
#Test the consistency by using a differenced χ2 consistency test,
#taking into account the uncertainties in the laboratory wavelengths of these lines
#to do the chisquared test we need the variance of g
##var_g= var_velocityabsone +var_velocityabstwo
#use propogation of errors
#variance of g
sig2_g= ((light/lambdarestone)**2)*unc[1][1] + ((light*param[1]/lambdarestone**2)**2)*(0.001)**2 +((light/lambdaresttwo)**2)*unc2[1][1] + ((light*param2[1]/lambdaresttwo**2)**2)*(0.001)**2
According to Morton 1991, the labatory wavelenths of these lines are 0.001.
Because I didn't want to latex my hand written work for propogation of errors, see canvas dropbox submission for written work for this part.
#using notes from lecture 18 on canvas
#finding the t statistic!
#assume mu is zero because doing difference in means test
tstat = (centvelone-centveltwo)**2/sig2_g
print(tstat)
from scipy.stats import chi2
#fidning PTE for t Test
#DOF is N-1, N= 2 because have two velocity points
DOF = 1
PTE = 1 - (chi2.cdf(tstat,DOF))
print(PTE)
10.059830398142388 0.001515371760480999
When testing the consistency of the two central velocities of the two absorption lines using a differenced $χ^2$ consistency test, taking into account the uncertainties in the laboratory wavelengths of these lines, I got the PTE=0.002 and $χ^2$ test statistic = 10.06. Based on lecture 18 notes, given that the PTE is low and the $χ^2$ test statistic is relatively high, we can assume that the two central velocities of the two absorption lines are inconsistent with each other or are different from each other.
Because I didn't want to latex my hand written work, see canvas dropbox submission for written work for this question.